#require(SLqPCR)
a <- read.csv("rat.csv",row.names=1,header=T)
sample <- as.factor(c(rep("Control", 3), rep("one", 3)))
res.Control <- selectHKgenes(a[sample=="Control",], method = "Vandesompele", geneSymbol = names(a), minNrHK = 2, trace = T, na.rm = T)
res.one <- selectHKgenes(a[sample == "one",], method = "Vandesompele", geneSymbol = names(a), minNrHK = 2, trace = TRUE, na.rm = TRUE)
ranks <- data.frame(c(1, 2:4), res.Control$ranking, res.one$ranking)
names(ranks) <- c("rank", "Control", "100 ppm")
print(ranks)
#write.csv(c(res.Control$mean,res.one$mean),file = "Gene stability.csv")

#install RColorBrewer from CRAN > require(RColorBrewer) 
#stability of gene expression upon exclusion of less stable genes
mypalette <- brewer.pal(4, "Set1")
matplot(cbind(res.Control$meanM, res.one$meanM), type = "b", ylab = "Average expression stability M", xlab = "Number of remaining control genes", axes = FALSE, pch = 19, col = mypalette, ylim = c(0, 0.025), lty = 1, lwd = 2, main = "Gene stability enamel organ")
axis(1, at = 1:3, labels = as.character(4:2))
axis(2, at = seq(0.000, 0.025, by = 0.01), labels = as.character(seq(0, 0.025, by = 0.01)))
box()
abline(h = seq(0, 0.025, by = 0.01), lty = 2, lwd = 1, col = "grey")
legend("topright", legend = c("No treatment","100 ppm"), fill = mypalette)

#variation of expression upon inclusion of additional genes beyond 3 in the calculation
mypalette <- brewer.pal(4, "PiYG")
barplot(cbind(res.Control$variation, res.one$variation), ylab="Pairwise variation",ylim=c(0,0.01),beside = TRUE, col = mypalette, space = c(0, 2), names.arg = c("No treatment", "100 ppm"),main="Variation of gene expression enamel organ")
legend("topright", c("V3/4", "V2/3","V3/4", "V2/3"), xpd = TRUE, horiz = F, inset = c(0, 0), fill = mypalette, ncol=2)
abline(h = seq(0, 0.01, by = 0.003), lty = 2, col = "grey")
#write.csv(c(res.Control$variation,res.one$variation),file = "Variation.csv")
